subroutine peo_init()
use molecule
implicit none
   !-------------------------------------------------------------------------
   character (len=4) :: HO_e(2), CH2CH2O_e(7), CH2CHFO_e(7), CH2CH3_e(7)
   real(q)           :: HO_p(3, 3), CH2CH2O_p(3, 8), CH2CHFO_p(3, 8), CH2CH3_p(3, 7)
   real(q)           :: HO_c(2), CH2CH2O_c(7), CH2CH2O_c_f(7), CH2CH2O_c_l(7), CH2CH2O_c_nf(7), CH2CHFO_c(7), CH2CH3_c(7)
   integer           :: HO_n(4, 2), CH2CH2O_n(4, 7), CH2CHFO_n(4, 7), CH2CH3_n(4, 7)
   integer           :: HO_rd(2), CH2CH2O_ld(2), CH2CH2O_rd(2), CH2CHFO_ld(2), CH2CHFO_rd(2), CH2CH3_ld(2)
   integer           :: nx, ny, nz, nn, id, i, j, ii, jj, href
   integer           :: dleft, dright, nleft, nright
   real(q)           :: hrefpos(3)
   logical           :: First_EO = .true., Follow_FEO = .false.
   character(len=100):: fmtstr
   !-------------------------------------------------------------------------
   HO_e  = (/ "HO", "OH" /)
   HO_c  = (/ 0.4110, -0.6844 /)

   CH2CH2O_e = (/ "C ", "H ", "H ", "C ", "H ", "H ", "O " /)
   CH2CH2O_c   = (/ 0.0023, 0.0935, 0.0935, 0.0023, 0.0935, 0.0935, -0.3786 /)
   CH2CH2O_c_f = (/ 0.1097, 0.1130, 0.1130,-0.0425, 0.0850, 0.0850, -0.3786 /) ! Charges for first -CH2CH2O-
   CH2CH2O_c_l = (/ 0.0018, 0.0931, 0.0931,-0.0012, 0.0902, 0.0902, -0.3786 /) ! Charges for last  -CH2CH2O-
   CH2CH2O_c_nf= (/-0.0049, 0.1060, 0.1060, 0.0050, 0.0935, 0.0935, -0.3786 /) ! Charges for -CH2CH2O- next to CH2CHFO-

   CH2CHFO_e = (/ "C ", "H ", "H ", "C ", "H ", "F ", "O " /)
   CH2CHFO_c   = (/-0.0221, 0.1060, 0.1060, 0.2253, 0.1160,-0.1731, -0.3786 /)

   CH2CH3_e  = (/ "C ", "H ", "H ", "C ", "H ", "H ", "H " /)
   CH2CH3_c  = (/ 0.0002, 0.0773, 0.0773, -0.2453, 0.0969, 0.0969, 0.0969 /)

   HO_n  = reshape( (/ 2, 0, 0, 0, 1, -1, 0, 0 /), (/4, 2/))
   HO_rd = (/ 2, 2 /)

   CH2CH2O_n = reshape( (/-1, 2, 3, 4, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, 7, 4, 0, 0, 0, 4, 0, 0, 0, 4, -1, 0, 0 /), (/4, 7/))
   CH2CH2O_ld = (/ 1, 1 /)
   CH2CH2O_rd = (/ 7, 2 /)

   CH2CHFO_n = reshape( (/-1, 2, 3, 4, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, 7, 4, 0, 0, 0, 4, 0, 0, 0, 4, -1, 0, 0 /), (/4, 7/))
   CH2CHFO_ld = (/ 1, 1 /)
   CH2CHFO_rd = (/ 7, 2 /)

   CH2CH3_n  = reshape( (/-1, 2, 3, 4, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, 7, 4, 0, 0, 0, 4, 0, 0, 0, 4, 0, 0, 0 /), (/4, 7 /))
   CH2CH3_ld = (/ 1, 1 /)

   ! Position of H, O, and expected bonding atom of O
   HO_p      = reshape( (/ -0.0994D0, -0.7841D0, -0.5447D0, 0.D0, 0.D0, 0.D0, 1.2954D0,  0.0000D0,  0.6027D0 /), (/3, 3 /))

   ! Position of CH2CH2O and expected bonding atoms to O
   CH2CH2O_p = reshape( (/  0.0000D0,  0.0000D0,   0.0000D0, &
                           -0.2285D0, -0.8951D0,   0.6124D0, &
                           -0.2286D0,  0.8944D0,   0.6128D0, &
                            1.4334D0,  0.0003D0,  -0.5404D0, &
                            1.6240D0, -0.8876D0,  -1.1726D0, &
                            1.6243D0,  0.8890D0,  -1.1721D0, &
                            2.4600D0,  0.0000D0,   0.5662D0, &
                            3.7712D0,  0.0000D0,  -0.0000D0 /), (/3, 8/))

   ! Position of CH2CHFO and expected bonding atoms to O
   CH2CHFO_p = reshape( (/  0.0000D0,  0.0000D0,   0.0000D0, &
                           -0.2285D0, -0.8951D0,   0.6124D0, &
                           -0.2286D0,  0.8944D0,   0.6128D0, &
                            1.4334D0,  0.0003D0,  -0.5404D0, &
                            1.6240D0, -0.8876D0,  -1.1726D0, &
                            1.9000D0,  0.9000D0,  -1.2500D0, &
                            2.4600D0,  0.0000D0,   0.5662D0, &
                            3.7712D0,  0.0000D0,  -0.0000D0 /), (/3, 8/))

   ! Position of CH2CH3
   CH2CH3_p  = reshape( (/  0.0000D0,  0.0000D0,   0.0000D0, &
                           -0.2285D0, -0.8951D0,   0.6124D0, &
                           -0.2286D0,  0.8944D0,   0.6128D0, &
                            1.4334D0,  0.0003D0,  -0.5404D0, &
                            1.6240D0, -0.8876D0,  -1.1726D0, &
                            1.6243D0,  0.8890D0,  -1.1721D0, &
                            1.7600D0,  0.0000D0,   0.5662D0 /), (/3, 7/))
   !-------------------------------------------------------------------------
   write(*, '(10x, "We are going to create a molecule of " )')
   write(*, '(10x, "       HO-[(CH2CH2O)x-(CH2CHFO)y]n-(CH2CH2O)z-CH2CH3" )')
   write(*, '(10x, "Please input the desired x, y, n, z (4 int): ", $)')
   read(*, *) nx, ny, nn, nz
   if (nx.lt.0) nx = 2
   if (ny.lt.0) ny = 0
   if (nz.lt.0) nz = 0
   if (nn.lt.0) nz = 1
   write(*, '(10x, "The desired molecule is: ", $)')
   fmtstr = '("HO-(CH2CH2O)", I??, $)'
   write(fmtstr(19:20), '(I2.2)') int(log10(dble(nx)+0.1))+1
   write(*, fmtstr) nx
   fmtstr = '("-(CH2CHFO)", I??, $)'
   write(fmtstr(17:18), '(I2.2)') int(log10(dble(ny)+0.1))+1
   if (ny.gt.0) write(*, fmtstr) ny
   fmtstr = '("-(CH2CH2O)", I??, $)'
   write(fmtstr(17:18), '(I2.2)') int(log10(dble(nz)+0.1))+1
   if (nz.gt.0) write(*, fmtstr) nz
   write(*, '("-CH2CH3")')
   if (((nx + nz).lt.2) .or. ((nx+nz).lt.3.and.ny.gt.0)) then
      write(*, '(/,10x, "Please pay attention to the charge neutrality of the molecule!")')
   endif
  
   natom    = 9 + 7*((nx + ny)*nn + nz)
   neighmax = 4
   molname  = "PEO"
   !
   if ( allocated(neigh   )) deallocate( neigh )
   if ( allocated(element )) deallocate( element )
   if ( allocated(molpos  )) deallocate( molpos  )
   if ( allocated(chg     )) deallocate( chg     )
   allocate( element(natom), neigh(neighmax, natom), molpos(3,natom), chg(natom) )
   !
   id = 0
   ! HO
   do ii = 1, 1
      do i = 1, 2
         id = id + 1
         element(id)   = HO_e(i)
         chg(id)       = HO_c(i)
         neigh(:, id)  = HO_n(:, i)
         molpos(:, id) = HO_p(:, i)
      enddo
      hrefpos = HO_p(:, i)
      dright = HO_rd(1)
      nright = HO_rd(2)
      !
   enddo
   !
   do jj = 1, nn
      ! 1st -CH2CH2O-
      do ii = 1, nx
         href  = id
         dleft = href + CH2CH2O_ld(1)
         nleft = CH2CH2O_ld(2)
         neigh(nright, dright) = dleft
         do i = 1, 7
            id = id + 1
            element(id) = CH2CH2O_e(i)
            if (First_EO) then
               chg(id)     = CH2CH2O_c_f(i)
            else if (ii == 1 .and. Follow_FEO) then
               chg(id)  = CH2CH2O_c_nf(i)
            else if (jj == nn .and. ii == nx .and. ny + nz < 1) then
               chg(id)     = CH2CH2O_c_l(i)
            else
               chg(id)     = CH2CH2O_c(i)
            endif
            where (CH2CH2O_n(:, i).gt.0) neigh(:, id) = href + CH2CH2O_n(:, i)
            where (CH2CH2O_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = hrefpos + CH2CH2O_p(:, i)
         enddo
         neigh(nleft, dleft)   = dright
         dright = href + CH2CH2O_rd(1)
         nright = CH2CH2O_rd(2)
         hrefpos = hrefpos + CH2CH2O_p(:, i)
         !
         First_EO = .false.
         Follow_FEO = .false.
      enddo
      !
      ! -CH2CHFO-
      do ii = 1, ny
         href  = id
         dleft = href + CH2CHFO_ld(1)
         nleft = CH2CHFO_ld(2)
         neigh(nright, dright) = dleft
         do i = 1, 7
            id = id + 1
            element(id) = CH2CHFO_e(i)
            chg(id)     = CH2CHFO_c(i)
            where (CH2CHFO_n(:, i).gt.0) neigh(:, id) = href + CH2CHFO_n(:, i)
            where (CH2CHFO_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = hrefpos + CH2CHFO_p(:, i)
         enddo
         neigh(nleft, dleft)   = dright
         dright = href + CH2CHFO_rd(1)
         nright = CH2CHFO_rd(2)
         hrefpos = hrefpos + CH2CHFO_p(:, i)
         !
         Follow_FEO = .true.
         !
      enddo
      !
   enddo ! Loop on -[(CH2CH2O)-(CH2CHFO)]-
   !
   ! 2nd -CH2CH2O-
   do ii = 1, nz
      href  = id
      dleft = href + CH2CH2O_ld(1)
      nleft = CH2CH2O_ld(2)
      neigh(nright, dright) = dleft
      do i = 1, 7
         id = id + 1
         element(id) = CH2CH2O_e(i)
         if (First_EO) then
            chg(id)  = CH2CH2O_c_f(i)
         else if (ii == 1 .and. Follow_FEO) then
            chg(id)  = CH2CH2O_c_nf(i)
         else if (ii == nz) then
            chg(id)  = CH2CH2O_c_l(i)
         else
            chg(id)  = CH2CH2O_c(i)
         endif
         where (CH2CH2O_n(:, i).gt.0) neigh(:, id) = href + CH2CH2O_n(:, i)
         where (CH2CH2O_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = hrefpos + CH2CH2O_p(:, i)
      enddo
      neigh(nleft, dleft)   = dright
      dright = href + CH2CH2O_rd(1)
      nright = CH2CH2O_rd(2)
      hrefpos = hrefpos + CH2CH2O_p(:, i)
      !
      First_EO = .false.
      Follow_FEO = .false.
   enddo
   !
   ! 
   do ii = 1, 1
      ! -CH2CH3
      href  = id
      dleft = href + CH2CH3_ld(1)
      nleft = CH2CH3_ld(2)
      neigh(nright, dright) = dleft
      !
      do i = 1, 7
         id = id + 1
         element(id) = CH2CH3_e(i)
         chg(id)     = CH2CH3_c(i)
         where (CH2CH3_n(:, i).gt.0) neigh(:, id) = href + CH2CH3_n(:, i)
         where (CH2CH3_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = hrefpos + CH2CH3_p(:, i)
      enddo
      neigh(nleft, dleft)   = dright
   enddo
   !
   !-------------------------------------------------------------------------
end subroutine peo_init
